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Population dynamics in systems composed of cyclically competing species has been of increasing 
interest recently. Here, we investigate a system with four or more species. Using mean field theory, 
we study in detail the trajectories in configuration space of the population fractions. We discover 
a variety of orbits, shaped like saddles, spirals, and straight lines. Many of their properties are 
found explicitly. Most remarkably, we identify a collective variable which evolves simply as an 
exponential: Q oc e^*, where A is a function of the reaction rates. It provides information on the 
state of the system for late times (as well as for t — >■ —oo). We discuss implications of these results 
for the evolution of a finite, stochastic system. A generalization to an arbitrary number of cyclically 
competing species yields valuable insights into universal properties of such systems. 
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INTRODUCTION 



In the study of population dynamics for biological sys- 
tems or chemical reactions, it is customary to start with 
a set of ordinary differential/difference equations which 
model the time evolution of populations or concentra- 
tions. These would be appropriate if the various con- 
stituents of the system are well mixed (so that spatial 
structures can be ignored) and if the stochastic aspects 
of the evolution can be overlooked. In this sense, the 
term 'mean field theory' (MFT) is most appropriate, as 
both spatial and temporal fluctuations are absent. De- 
spite these limitations, MFT has proven valuable in sev- 
eral respects. As in other areas of statistical physics, 
it typically provides an adequate description in much of 
the space of (control) parameters, e.g.. Landau's theory 
of phase transitions. In addition, it yields significant in- 
sights into many interesting phenomena associated with 
the underlying non-linear dynamics. A well known exam- 
ple is the simple logistic map, which gives rise to a rich 
structure of bifurcation, chaos, and universality [H-Q. Fi- 
nally, it provides a stage on which hidden symmetries can 
be showcased readily. 

In this context, we are motivated to consider MFT 
for a simple model of population dynamics, involving 
four species competing cyclically. This is a generalization 
of systems with three cychcally competing species, also 
known as the rock-paper-scissors games ,which attracted 
considerable attention in recent years |8l-[29|. Labeling 
the species a, b, c, d, we allow a to 'prey on' b with some 
rate, and 6, c, d to 'prey on' c, d, a, respectively, with some 
other rates. Unlike earlier studies, we will not restrict 
our rates, but consider them in general. As we will show, 
apart from having a much larger and complex parameter 
space, the four species system (4SS) is qualitatively dif- 
ferent. Indeed, in MFT, the evolution of the three species 
system (3SS) is relatively trivial, since a non-linear in- 
variant (a hidden symmetry) renders all trajectories into 



closed orbits on a plane. As a result, there is little hint of 
the presence of absorbing states, which are of fundamen- 
tal interest for the study of extinction probabilities in a 
stochastic system. By contrast, such a scenario is not the 
case, typically, for the 4SS. In addition, when one species 
becomes extinct, our problem reduces to an apparently 
trivial limit of the 3SS. Yet, the MFT offers useful and 
quantitative predictions. In a separate publication 30], 
we will present extensive simulations of the full stochastic 
model and show that many of the MFT predictions are 
well born out. In this article, we will focus only on the 
deterministic MFT and the rich variety of phenomena it 
predicts. Though some of the preliminary results on the 
4SS have been presented earlier [3l| , this paper will delve 
into not only details of the analysis, but also provide new 
insights into the behavior of systems with S > A species. 

Before we begin, we should mention that there are 
studies for the general S case, though far few in number 
than for the 3SS. Most of these studies considered sys- 
tems on one- or two-dimensional lattices 
Due to the presence of spatial structure, these systems 
display far more complex behavior and so, the authors 
typically restrict themselves to special sets of rates (e.g., 
all equal, or all but one being equal). Finally, some stud- 
ies, motivated by real-world systems, focused on a large 
number of competing species with more complicated in- 
teractions schemes |39|, 140]. Here, we will focus on phe- 
nomena that emerge from arbitrary (positive) rates, but 
consider only systems with no spatial structure. An inter- 
esting application of a four-species model without spatial 
structure to the endogenous and exogenous origins of dis- 
eases can be found in (4]| . 

This paper is organized as follows. Though our main 
focus will be MFT, we will first discuss how the differen- 
tial equations arise from a fully stochastic, microscopic 
model of the 4SS. Details of the latter (which can be 
simulated by Monte Carlo methods) and the associated 
master equation will be presented in Section II. The ap- 
proximation scheme which leads to the MFT will also 
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be provided. Section III will be devoted to the analy- 
sis of the MFT equations, with a number of explicit re- 
sults. Their implications for the stochastic evolution will 
be considered in Section IV. By generalizing to systems 
with 5* > 4 cyclically competing species, we gain deeper 
insights into how certain aspects of the solutions are 'uni- 
versal.' These considerations will be presented in Section 
V. Readers who are comfortable with abstract formula- 
tions will find that some of the conclusions in Sections 
III are just special cases. Finally, we will summarize our 
findings in a concluding section and sketch possible av- 
enues for future research. 

Interestingly, we found [3l[ that, unlike in the 3SS, the 
four species form partner-pairs, much like in the game 
of Bridge. In a system with N individuals, there are 
2 {N + 1) absorbing states, most of which consist of a 
surviving partner-pair. Since MFT is best suited for sys- 
tems with large N, we should expect that all the final 
states in this approximation to consists of a coexisting 
pair: either a-c or b-d. Of course, many interesting issues 
associated with the full, stochastic model (e.g., extinc- 
tion probabilities) cannot be answered by MFT, while 
simulations reveal complex extinction scenarios that de- 
pend on both the predation rates and initial conditions. 
Since those investigations are quite extensive, they will 
be reported elsewhere [30| . 



have 2 (A'' -I- 1) absorbing states here instead of just 3 (re- 
gardless of N). We should also remark that each face of 
the tetrahedron is also 'absorbing' in the sense that tran- 
sitions into the face are irreversible. Within each face, 
the problem is a special limit of the 3SS, namely, one of 
the three rates being zero. 

Given the interaction rules, the master equation for 
P{{Nm}',T), the probability for finding the system t 
steps after a particular initial distribution Pq {{Nm}), 
can be easily written. It provides the change in P over 
one step: 

Pi{N^};T+l)-Pi{N„,};T) 

Pa {Na - 1) im + 1) 



N{N ~1) /2 

Pb (Nb - 1) (iVe + 1) 

N{N-l)/2 

p, {N, - 1) {Ng + 1) 

N{N-l)/2 
Pd {Nd - 1) {Na + 1) 



N{N -1) /2 
^ -P({iV„};r) 



N[N -I) /2 



P{Na-l,Nk + l,N,,Nd]T) 
P{Na,Nb-l,N,+ \,Nd\T) 

P{Na,Nb,Nc-l,Nd + l;T) 
P{Na + l,Nb,N,,Nd~l;T) 

(1) 



where 



II. THE MICROSCOPIC STOCHASTIC MODEL 
AND ITS MEAN FIELD APPROXIMATION 

Our four species system consists of 7V„i individuals of 
species m = (a,6, c, d). With no spatial structure, any 
individual can interact with any other, in a cyclically 
competing manner. In a time step of a fully stochastic 
model, a random pair is chosen and, if they are ' cyclically 
different,' the following interactions occur with probabil- 
ities pm ■ 

a + b ^ a + a; b + c ^ b + b 
c + d ^ c + c; d + a ^ d + d 

Thus, the total number in the system 

m 

is a constant and the full configuration space is just a set 
of points within a regular tetrahedron (of length N on 
each side, with the four vertices being N = Nm for one 
of the rn's). Of course, we will later consider the fractions 
Nm/N, which are natural variables in the large N limit 
and the mean field approximation. 

Let us emphasize that ac and bd pairs do not interact, 
so that it is possible for the final (absorbing) state to dis- 
play coexistence of these pairs. In other words, any point 
along the a-c and b-d edges of the tetrahedron represents 
such a state. Indeed, this property provides the first ma- 
jor difference between our 4SS and the simpler 3SS: We 



Z = PaNaNb + PbNbNc + PcNcNd + PdNdNa . (2) 

For the initial Pq ({-^m}), it is sufficient to use 5 distribu- 
tions. Indeed, every simulation run begins with a single 
point (e.g., the symmetry point Nm = N/A)\ Since the 
master equation is linear, the evolution starting with any 
other distribution is just the sum of the P's that begin 
as 5's at various appropriate points. 

From the master equation for P {{Nm.} ] t), we can eas- 
ily derive a partial differential equation for the associated 
generating function. But, to find the solution for either 
equation is far from trivial. Instead, we turn to the mean 
field approximation, which should be adequate for large 
{Nm}- Thus, we can expect break downs near the ex- 
tinction of one or more species, i.e., near an 'absorbing 
face' or an absorbing state. 

In a MFT, the focus is the evolution of the averages of 
say, the fractions Nm/N, denoted by 

A{r) EE (Na/N)^^ {Na/N)P{{Nm};r) , 

B{t) EE {Nb/N).^,C{T)^{NjN)^,D{T)^{Nd/N)^ 

Of course, we have the constraint 

A + B + C + D = 1 (3) 

so that there are only three independent variables. 

Following standard routes, we can derive an equation 
for the changes, A{t + 1) — A(r), etc., from the mas- 
ter equation ([I} . These will involve averages of products 
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(e.g., (NaNi,)^) on the right. The mean field approxima- 
tion consists of neglecting all correlations and replacing 
the averages of products by the products of averages, so 
that the end result is a closed set of deterministic equa- 
tions for A{t), B (t), etc. Finally, by taking the N oo 
limit, defining t = t/N (which can be regarded as con- 
tinuous), and letting A{t + 1) ~ A{t) (1/iV) dtA (i), 
we arrive at the MFT equations: 



dtA 
dtB 
dtC 
dtD 



[kaB - kdD] A 
[kbC - kaA] B 

%D - kbB] C 
[kdA - kcC] D 



(4) 
(5) 
(6) 
(7) 



In these time units, the 'rates' km can be related to the 
original probabilities via km = 2pm- In the literature, 
however, time is often rescaled again, so as to normalize 
the rates to 



ka + kb + kc + kd = I 



(8) 



as well. 



III. ANALYSIS OF MFT EQUATIONS 

This section will be devoted to the implications of eqns. 
(jUl?]), assuming the initial fractions are {Aq, Bo,Co, Dq). 
Given the normalization ([5]), our parameter space is, in 
general, quite large: 6-dimensional. We will begin with a 
restricted problem - evolution on one of the faces of the 
tetrahedron - which is explicitly solvable. 



Meanwhile, B ^ I - A - C and A^''" (C^'^") is given 
by RV {R/V), so that the full solution for V {t) can be 
obtained through the integral 



V{t) 



dv 



Vo 



V 1 - [Rvf'^''' - {R/vf'^''- 



2kakbt (12) 



where Vq = Ag^'Cj^'^". Though analytically precise, this 
form of the solution does not provide much insight on 
the evolution of the system. Instead, a schematic plot 
of sample trajectories, as shown in Fig. [l] should be 
much more helpful. Representing various values of R, 
the curved lines (dashed, red online) are generalized hy- 
perbolas ([9]). The initial fractions {Aq, Bo,Co) dictate 
which curve that system will follow, while the evolution 
takes it towards larger A (in the direction of the arrows, 
blue online), ending at a point {Af,0,Cf) on the a-c line. 
Note that the straight line (dot-dashed, green online), 
given by kbC = kaA, intersects each of these curves at 
right angles, corresponding to B being maximal at these 
points. Finally, by exploiting the invariant i?, the values 
Af and C/ (= 1 — Af) can be found from a (generally 
transcendental) equation, e.g., 

A''/ (l-Aff'' = A[;^Co^" . 

Since ka,kb G (0,1), a schematic plot of the left hand 
side reveals that there are two solutions. Of these, the 
larger one is the final Af, since dtA > 0. Of course, 
numerical methods will be needed typically for obtaining 
the explicit Af. These considerations are useful when 
MFT predictions are compared to simulation data for 
times after one of the four species goes extinct [sO, HH . 



A. Evolution when one species is extinct 

Clearly, this restricted problem is not only much sim- 
pler than the full 4SS, it is also considerably simpler than 
the general 3SS. The reason is obvious: One of the three 
remaining species does not 'consume' anyone and so, its 
numbers can never increase. Without loss of generality, 
let us consider the face with D = with arbitrary initial 
{Aq, Bq, Co). In this case, neither kc nor kd plays a role, 
while eqns. (I4][7l) reduce significantly. As a result. 



(9) 



is an invariant (similar to that in the full 3SS [20|), so that 
it is always fixed at Aq''Cq'^ . Of course, our problem now 
reduces to one with a single variable, which we choose to 
be 



V = A^^C-^-^ 



(10) 



so that its evolution is monotonically increasing until B 
vanishes: 



dtlnV ^ 2kakbB > 



(11) 



B. A general criterion for survival/extinction 

Turning to the full 4SS, we identify a single parameter 
that determines which of the pairs survive in the long 
time limit. To find this criterion, we note that eqns. (j4]- 
[7]) can be written as 

dtlnA = kaB - kdD; dtlnC = k^D - kbB (13) 
dtlnB = kbC-kaA; dtlnD ^ kdA ~ k^C (14) 

since none of the fractions vanishes in finite t. This form 
is quite natural, since exponential growth/decay in pop- 
ulations is common. It also exposes clearly the coupling 
between the pairs a-c and b-d. Constructing appropriate 
linear combinations, we showcase the contributions from 
a single species to the growth/decay of two other species: 

dt[kblTiA + kalnC] = XD 

dt [fcclnA + fcrflnC] = \B 

dt[kchiB + kb\nD] = -\A 

dt[kdlTLB + ka\TLD] = -AC 
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FIG. 1: (Color online) Schematic plot of trajectories in the 
{A,B,C) plane in case species d dies first. The dashed red 
lines show the generalized hyperbolas ((Ojl along which the sys- 
tem evolves by increasing A. At the intersections of the hy- 
perbolas and the dot-dashed green lines B is maximal. 

where we have introduced the key control parameter 

A EE kakc ~ hkd ■ (15) 

Adding and subtracting appropriately, we see that the 
quantity 

'5 = Bka+kd Dka+kt ^^^^ 

evolves in an extremely simple manner: 

Q{t) = Q{0)e^' . (17) 

The form of Q also points to a simple interpretation of 
the behavior of the species. As in the card game, Bridge, 
a-c and b-d form opposing partner-pairs. For example, a's 
action (consuming b) benefits c, while c's action benefits 
a. Thus, we may refer to a-c and b-d as partners. Since no 
fraction can exceed unity, the only way for Q to vanish 
or diverge is for one of the pairs to go extinct. Thus, 
the sign of A is the selection criterion for, and Q is a 
quantitative measure of, which pair survives. 

In this connection, we see that, unlike the 3SS [l^l, the 
'weakest' is not necessarily the survivor (or among the 
survivors) here. Indeed, if both members of a partner 
pair are weak (compared to their opponents), then the 
sign of A predicts the eventual demise of this pair, while 
the magnitude of A hints at how rapidly they will die out. 
As we have pointed out 31], the general maxim seems to 
be: "The prey of the prey of the weakest is the least likely 
to survive." In the 3SS, this maxim also applies and is 
consistent with 'survival of the weakest.' There, 'the prey 
of the prey of the weakest' happens to be its predator, 
the demise of which naturally enhances its survival. 



C. Saddle shaped orbits and fixed points for 
systems with kakc = ki,kd 

Clearly, special properties will appear in a 4SS with 
A = 0, when Q becomes an invariant. Indeed, there are 
two invariants [1, [2^ [3l[ . Neither pair goes extinct and 
the system evolves along periodic, closed orbits in con- 
figuration space. It is tempting to regard Q as the gen- 
eralization of the quantity R = A*'* B^''C'^'^ (introduced 
in (20| ) in the 3SS, except that R is invariant for any 
set of rates. As will be shown in Section V, there are 
also fundamental differences between systems with even 
and odd number of cyclically competing species. As a re- 
sult, direct comparisons between R and Q are not helpful, 
though both play important roles in identifying collective 
degrees of freedom that evolve trivially. 

Proceeding to investigate the extra invariant besides 
Q, we find an appealing, symmetric way to display the 
two constants of motion: 

f^Akbfjk^ and g = B^''D^-^ (18) 

(or equivalently, /' = A^-'C^'^ and g' = B^'=D^>'). 
These generalize the products AC and BD in [6] , where 
the rates are all unity. Since invariants are fixed by the 
initial conditions (/ = Aq^'Cq", etc.), it is natural to de- 
fine the variables 

^ {A/A^f-; p^^iC/Co)'- (19) 
p, ^ (B/Bo)"' ; Pd ^ (D/Do)"^ (20) 

Then, the constants of motion are simply the products: 
PaPc = 1; PbPd = 1 • (21) 

They define hyperbolic sheets through the tetrahedron 
and their intersection is a closed loop that resembles (the 
rim of) a saddle. Fig. [5^ shows an example of such an 
orbit, for the case ka = 0.4, ki, — 0.4, kc = 0.1, kd — 
0.1. Here, the MFT eqns. (|4][7]) have been integrated 
numerically using a standard fourth-order Runge-Kutta 
scheme. Fig. shows the ever-lasting oscillations in 
A, ...,D associated with this orbit. 

Many characteristics of a closed orbit in our tetrahe- 
dron are accessible analytically. We will summarize the 
findings here and defer the details to the Appendix. Let 
us start with the simplest: the average position (over one 
period, T): 

o I r ° I F 

A=- A {t) dt- B^- B {t) dt; etc. 
^ Jo ^ Ja 

Integrating eqns. (|13I14|) . we find that the left vanish 
identically. Thus, we conclude 

ka A= kb C; ka B^kdD ■ (22) 

Below, we will discuss the presence of a line of fixed- 
points for A = systems, its properties given by the rates 
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FIG. 2: (Color online) (a) Example of a closed loop (solid 
curve) in the tetrahedron forming the configuration space, 
encircling the line of fixed points (green dashed line). Note 
that the fixed line bridges absorbing states on the a-c and 
b-d edges. The loop is generated by solving the mean field 
equations with a fourth order Runge-Kutta scheme, using 
time step t = 0.00001. The rates (fca, fci, fee, fed) used are 
(0.4,0.4,0.1,0.1), yielding A = 0. (b) Evolution of the four 
species fractions as a function of time. The initial fractions 
{Ao, Bo, Co, Do) are (0.02, 0.10, 0.48, 0.40). Next to each curve 
is its species label. The red and blue open circles mark the 
values of B and D, respectively, associated with the two turn- 
ing points of A, located by the black arrows. Note that both 
red (blue) circled values are the same. See text for details. 



solely. It is perhaps not surprising that yA, B, C, Dj lies 

on this line. Its precise location depends on the initial 
conditions, but finding this explicit dependence remains 
elusive. Another conclusion from this consideration is 
that all closed orbits enclose this non-trivial fixed line. 

A more detailed description of our closed orbit is its 
projection onto, say, the pa-pd plane. One possible rep- 
resentation is to exploit eqns. (jl9l20l3|) : 



AaPa 



1 



(23) 



In the Appendix, we will show how this ungainly looking 
equation can be transformed into the familiar equation 
for a circle: + 0^ = 1. We can gain some insight into 
this orbit by studying the extremal points. In particular, 
focusing on A and D of an orbit, we denote the extremes 
through the expressions: 



These extremal points are given by the solution to the 
following equations: 



A± + JaA^ 



1-Ka{ B^'D^' 



l/(fca + fed) 



l/ika+kt) 



where J and K are constants: 



Ja = C'o^o''^'^"' Jd — BqDq° 

kd/{ka+kd) /, \kc,/(ka+kd) 

ka/(ka+kb) 



•\ka/kd 



kb/{ka+kb) 



(24) 
(25) 

(26) 
(27) 

(28) 



Being transcendental in general, these can be solved only 
numerically. Nevertheless, we can appreciate that there 
are typically two solutions to each (by a simple sketch of, 
say, x+x~P). Obviously, the smaller (larger) of the two is 
the minimum (maximum). Interestingly, when a species 
is at its extremum, the fractions of its opposing partner- 
pair assume the same value. For example, when A is at 
either turning point (^±), B takes on the same value. 



k'^-k-'^-B^JD. 



l/(ka + kd) 



value here: 



Similarly, D takes on one 

l/ika + kd) 



To illustrate, we 



-kd ofcd n^i^ 

^0 ^0 

mark one set in Fig. [^Jj with open circles. As will be 
discussed in Section IV, the minima {A_, B^,C^, D_) 
can play a role in predicting survivability in stochastic 
systems. 

Finally, like the 3SS, there are non-trivial fixed points. 
Indeed, with one more degree of freedom than 3SS, we 
find a line of such points here. Note that these are neu- 
trally stable, unlike the trivially time-independent, ab- 
sorbing states. The simplest route to this line is to set 
the right of eqns. (I3][71) to zero. Denoting the values 
on this line by {A* , B* ,C* , D*), we see that they sat- 
isfy kaA* = kbC* and kaB* = kdD* . In other words, 
the line is the intersection of the two planes. Exploiting 
A + B + C + D = 1, we find an elegant way to display 
these values. Introducing a parameter 7 S [0, 1], we have 



{A*,C*) 



(kb, ka) 
ka ~l" ki) 



r, {B\D* 



ka + kd 



(1-7) 



(29) 

Clearly, this line is straight, and bridges the two ab- 
sorbing states on the a-c and b-d lines, as illustrated in 
Fig. [5^. For the neighborhood of this line, it is natu- 
ral to study linearized versions of eqns. (|4][7l). Writing 
A = A*-\-Aa with Aa ■^A*,B = B*-\-Ab with At < B*, 
etc., we first note that A„ = and dt (A^ -I- Aj,) = 
= dt {Ab -I- Ad). Then, by defining q = kaAa - k^A^ 
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(«) (b) 



FIG. 3; (Color online) (a) A-C and (b) B-D plane. When a 
trajectory enters the unshaded region in (a) and the doubly 
shaded region in (b), both B and D decrease whereas at the 
same time both A and C increase. See main text. 



and p = kaAb — fc^A^, we can cast the linearized equa- 
tions in the form of a standard simple harmonic os- 
cillator: q = {ka + h)A*p; p = ~ {ka + k^) B*q. In 
other words, the saddle-like orbits become planar while 
th e system evolv es along ellipses, with angular frequency 
yjkbkd'y (1 — 7). Remarkably, contrary to natural expec- 
tations, the fixed line is not normal to the orbit plane in 
general. 

To summarize, if A 0, the evolution is relatively 
simple. Every starting point is associated with a closed, 
saddle-shaped orbit, on which the system remains for- 
ever. In general, analytic closed form expressions for 
them are not known (as far as we are aware). However, 
all enclose a fixed line and many of their properties can 
be computed analytically. Obviously, for finite systems 
evolving stochastically, these conclusions are valid only 
approximately - as long as the system is far from ab- 
sorbing states. 



D. Systems with kakc 7^ ki,kd • spirals and arrows 

With such rates, Q grows/decays exponentially, so that 
non-trivial fixed points cannot exist. Since all fractions 
are bounded by unity, this implies that either BD or AC 
vanishes in the large t limit. In other words, for a system 
with finite N , extinction of one of the species must occur 
quite rapidly, especially if A is large. Before we continue 
to the analysis of these cases, let us provide a simple and 
intuitive picture, for appreciating the significance of A. 

From eqns. (|13I14[) . we see that the rise and fall of 
each species are associated with the relative values of the 
opposing pairs. For example, we can display a line in 
the B-D plane {kaB = kdD, dashed line denoted by d 
in Fig. ^jp, on one side of which A increases (shaded re- 
gion, yellow online) and on the other, A decreases. Of 
course, this 'line' is actually a plane in our configuration 
space tetrahedron. When the orbit crosses this plane, A 
reaches a (local) extremum and turns around. A simi- 



A 




FIG. 4: (Color online) Typical orbit for A > that starts at 
the solid green circle and spirals towards an absorbing state, 
indicated by the red cross on the a-c edge. The data shown 
here have been obtained for rates ka ~ 0.45, kb = 0.33, fee = 
0.14, kd = 0.08 and initial fractions A = 0.02, B = 0.10, 
C = 0.48, D = 0.40. 



lar line {kbB = kcD; c in Fig. [Sjj) can be plotted for 
the 'divide' between increasing (shaded region, green on- 
line) and decreasing C. The other set of lines, b and d 
in the A-C plane, are shown in Fig. The significance 
of A is now clear. For A = 0, a = c, so that there are 
no regions where both A and C increase or decrease to- 
gether. Meanwhile, we also have b = d so that B and 
D also share such a property. By contrast, if A > 0, a 
'gap' opens between these pairs of lines in such a way 
that in the region between d and c, both A and C in- 
crease (e.g., overlap of shaded regions in Fig. [3}d, a solid 
'wedge' in the tetrahedron). On the other hand, in the 
region between b and d, both B and D decrease (un- 
shaded region in Fig. [3^). As a result, when the system 
gets into this domain (intersection of doubly-shaded and 
unshaded regions within our tetrahedron - i.e., an 'irreg- 
ular tetrahedron'), B and D will monotonically decrease, 
exponentially at late times. Simultaneously, A and C will 
monotonically increase, towards a final point {Af,Cf) on 
the a-c edge (solid diagonal line in Fig. |3^, red on line). 
In particular, the system will end on the heavy solid seg- 
ment (red on line) of this line. In Fig. 21 we illustrate 
these features with such a (A > 0) case, showing a typical 
open orbit (starting at the solid circle, green online), spi- 
ralling towards an absorbing state on the a-c edge (cross, 
red online). 

To appreciate such monotonic behavior better, we 
present here a unique trajectory. Instead of being a spi- 
ral, it is as 'straight as an arrow.' A system started on 
this (straight) line will move along it for all t. In fact, as 
A 0, this line becomes the fixed line above. Here 
we simply quote the result. Readers interested in the 
origins and general characteristics of such 'arrows' will 
find details in Section Vb, which is devoted to systems 
with arbitrary (even) S. We start with an Ansatz for the 
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fractions 



Ah {t), B[l-h (t)] ,Ch{t), D[l-h {t)] 



(30) 



where A, B, C, D are constants. Clearly, these points lie 
on a straight line in the tetrahedron. For h G [0, 1], this 
line joins the points (A, 0,(7, 0) and (0,5,0,1?) on the 
a-c and b-d edges, respectively. Inserting these into eqns. 
(|4|5p . we find that they can be solved provided 

dtlnh — uj {I — h) ; dt In [\ ~ h) — —ujh 

and 

cj = kaB - k^D = kcD - k^B = kaA - k^C = kcC - kdA 



Thus, the result is 
{A,B,C,D) = 

so that, explicitly. 



(ki) kc, kc -\- kd, ka '\- kfi, ka -i^ kh) 



ka ^ ki, -\- kc -\- kd 



(31) 



_ ka {h + kc) - kb {ka + kd) 
^ — ; '. — ; '. — ; — '. — ; l^o^ ) 



ka ~\~ kjj -\- kc -\- kd 
A 



fca + fcfe + fcc + kd 



(33) 



Since h satisfies dh/dt = ujh{l — h), we may choose 
h{0) = 1/2 and write the evolution along this straight 
trajectory explicitly: 



hit) 



1 



(34) 



Note that, indeed, h e [0,1] for all t. Further, as 
A — >■ 0, so does oj and h become 'frozen,' playing the 
role of 7 in ([29]) . Of course, in this limit, we easily see 
that, e.g., the ratio A/C = {kf, + kc) / [ka + kd) becomes 
h (1 + kd/ka) I (ka + kd) which is Ujka = A*/C*. It is 
reassuring that this 'arrow' trajectory for A ^ is con- 
tinuously linked to the fixed line in A = 0. Apart from 
this special line, all orbits (which we found numerically) 
spiral like a corkscrew, to various degrees, between the 
a-c and h-d edges. Classifying the characteristics of these 
spirals, a task well suited for future studies, will be highly 
non-trivial. 

Whether we consider the arrow or spiraling trajecto- 
ries, we remark that the time-reversed evolution (starting 
from any initial point) is also of interest. In this case, the 
opposite will occur: For systems with A > 0, say, Q will 
decrease monotonically and the system 'ends up' some- 
where on the opposite (h-d) edge. In this sense, the full 
dynamics for A 7^ can be regarded as a family of spirals 
around an arrow, bridging the two opposite edges The 
initial condition simply serves to pick out the appropri- 
ate orbit, while the evolution consists of following it to 
an absorbing edge. Shown in Fig. [S]is such a spiral, but 
with two branches (black line vs. magenta [grey] line) 
corresponding to < — > ±cx) (starting at the same point). 
Underlying this dynamics appears to be a kind of 'du- 
ality' symmetry (associated with A,t <^ — A, — t). Ex- 
ploring this aspect is likely a worthy pursuit, but clearly 
beyond the scope of this paper. 



D 



B 




C C 



(b) 



FIG. 5: (Color online) A typical 'arrow' and spiral, with A = 
-0.0273, from rates (ka, kt, kc, kd) = (0.35,0.42,0.09,0.14). 

(a) Starting at the symmetry point (all fractions being 0.25), 
the forward spiral is shown in magenta [grey]. The t — >■ — 00 
branch is shown as a black spiral. If started anywhere on the 
green arrow, the system will follow the arrow to the b-d edge. 

(b) A different perspective of (both branches of) the same 
spiral. 



IV. IMPLICATIONS FOR STOCHASTIC 
EVOLUTION 



When compared to the evolution of a fully stochastic 
4SS with finite N, the predictions of the deterministic 
MFT can be a good approximation for the average be- 
havior for cases with large iV^'s. Of course, MFT has se- 
rious limitations as well, apart from its inherent inability 
to describe fluctuations. The most glaring discrepancy 
lies with issues concerning extinction. Examples which 
readily come to mind are: What is the eventual extinc- 
tion probability of a species? At what time is a species 
most likely to go extinct? And what is the average time 
for extinction? Since the fractions in MFT never vanish 
in finite time, none of these can be answered. Never- 
theless, MFT can provide valuable insights [1^, some of 
which will be provided here. Extensive studies of the 
stochastic system and comparisons with MFT are be- 
yond the scope of this article and will be presented in a 
subsequent publication (soj . 

Let us consider first the A = case, where only closed 
orbits prevail and all species survive in the MFT. This be- 
havior resembles that in a 3SS. As pointed out in [20], the 
closest approaches a closed loop makes with each edge (of 
the configuration space triangle) are good indicators of 
the survival probabilities. We should caution the reader 
that, in the 3SS triangle, an edge is associated with the 
survival of a single species, as well as the extinction of 
one other species. Thus, it is natural to associate a min- 
imum distance to an edge (A™, in the notation of p(l ]) 
with the survival of species m (instead of the extinction 
of species m — 1). In our case, the extinction of a species 
induces the survival of the opposing partner-pair. Thus, 
we find it more convenient to use the language of extinc- 
tion of a species rather than likely survival of a pair. In 
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other words, the closest approach to one of the four faces 
(in our tetrahedron) is a good indicator of the extinc- 
tion probabihty of the associated species. For example, 
is the closest a loop comes within the face A = 0, 
which corresponds to a being extinct. Therefore, we can 
expect a species associated with the smallest in the set 
{A-, B-,C-, D-) to vanish first. In contrast to the 3SS 
case, the rates alone do not provide a useful criterion for 
extinction in our case, as we will show next. 

In the 3SS, the initial conditions identifies a closed or- 
bit and is associated with the invariant R = A*''' B^'^C^'^ . 
It is easy to find the closest approaches {A- , i?_ , C_ ) 
and to compute their ratios. For orbits with R 1, 
these ratios are well approximated by 

cx i?i/fc.-i/fcc ^ (35) 

B_/C_ cx i?i/fcc-i/fc<. ^ (3g) 

C_/A_ cx i?i/fca-i/fc. . (37) 

Notably, the proportionality coefficients depend only on 
the rates and not on R. Therefore, given a set of rates 
and a large population {N ^ 1), it is possible to ar- 
rive, through stochastic effects, at orbits where one of 
the fractions Nm/N — >■ <ti 1. Considering orbits 

with smaller and smaller i?, we see that the above ra- 
tios will vanish or diverge, depending on the rates alone. 
For example, if ka < kh,kc , then C_ will be the small- 
est as i? — )• 0, so that c is be the most likely to go 
extinct first. This argument leads to the conclusion in 
[2Qj: Extinction/survival probabilities are either zero or 
unity, in the limit of large N. By contrast, two invari- 
ants (/ = A''>'C^'',g = i3'=<i£)'=«) are associated with the 
closed loops in the 4SS. As a result, similar ratios of the 
closest approaches are not simple functions as in ([55]). 
Let us examine this issue further, by studying an ex- 
ample. We will again consider cases with small closest 
approaches, so that we can ignore the first term in eqn. 
(1^^ for finding, say, The results are 







1 




fed)" 
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+ kb) 


- -ka/ka 




^ 5V.. 


1 


-Kof 




- -ka/ka 



where Ka,d are given by (|27I28I) and are independent 
of /, 5. Even the ratios of each partner pair are not so 
clear-cut as in psp . since the 'other invariant' enters in 
a nontrivial way. For example, g appears in the ratio 

^ ^ fi/k,-i/ka _ -'^'"/'^^ _ 

Needless to say, 'cross-ratios' such as A^/B_ are even 
more intractable. In a stochastic evolution, / and g pre- 
sumably wander separately, so that it is a priori unclear 



which pair will win. Since the qualitative picture is al- 
ready quite complex, we will not pursue the route in [2^ 
where more quantitative arguments concerning scaling 
can be fruitful. Instead, devoting efforts towards the 
study of the full stochastics, along the lines of is 
likely to be more worthwhile. 

Turning to systems with A 7^ 0, we can reach simi- 
lar qualitative conclusions concerning how well the MET 
predicts the behavior of the stochastic model. Here, Q 
grows/decays exponentially and which pair survives is de- 
termined only by the sign of A. Phrased in stark terms, 
this conclusion states that neither its magnitude nor the 
initial configuration (as long as they lie within the tetra- 
hedron) plays a role. Clearly, this cannot be true for 
stochastically evolving, finite systems. From our sim- 
ulations (soj . we glean another maxim: "The prey of 
the prey of the strongest is quite likely to survive." Al- 
though this principle appears to be the same as the one 
posed above (albeit a double- negative version), the con- 
sequences for the end state can be drastically different. 
In particular, both |A| and the initial conditions, as well 
as the total population size, can reverse the predictions of 
MET. Let us illustrate with the example in [sij with 'ex- 
treme' rates: k„, = (0.1,0.001,0.1,0.7999). Since A > 0, 
a-c is predicted to win. When the initial population is 
set at Nm = (100, 700, 100, 100), 90% of the runs indeed 
end with the a-c pair. However, for a symmetric starting 
point of Nm = (100, 100, 100, 100), 97% of the runs end 
with b as the sole survivor, despite the initial Nb being 
much smaller! To accommodate these two very different 
outcomes, we apply the two different maxims. As b is 
the weakest and d is the strongest, each is the prey of 
the other's prey. In one case, the first maxim lead us to 
D vanishing first, leaving us with a-c coexistence. In the 
second scenario, d consumes a so rapidly that A vanishes 
first and b can only increase. But b is slow to consume 
c, so that d goes extinct next. As a result, b is the only 
survivor (rather than one of the more numerous states 
of b-d coexistence). To turn these qualitative considera- 
tions into quantitative conclusions, concerning the detail 
interplay between the rates and the initial configuration, 
will be a serious challenge. 

To summarize this section, let us reiterate that, while 
MET cannot predict all eventualities of a finite stochas- 
tic system, it can provide some hints. Specifically, the 
proximity of an MET trajectory to an absorbing face 
should be a good guide for the extinction probabilities 
of the associated species. For the 3SS, this guide led to 
very successful predictions, i.e., survival of the weakest. 
This success can be traced to the presence of only a few 
degrees of freedom, so that, along with the existence of 
the invariant R, the rates alone determine the relative 
proximity of the orbits to the three extinction lines. By 
contrast, the 4SS is sufficiently more complex that similar 
conclusions cannot be drawn. 
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SYSTEMS WITH S SPECIES 



A. Systems with odd S 



In this section, we present extensions of our findings 
to a system of N individuals consisting of an arbitrary 
number, S, of cyclically competing species. This gener- 
alization provides some insight into the general structure 
of the dynamics in such systems. In particular, these 
considerations lead us to the deeper origins of the ex- 
istence of collective variables like R and Q, as well as 
fixed points and 'arrows.' For any odd S, we show that 
there is precisely one (non-trivial) fixed point and one 
i?-like invariant. By contrast, for even S, a Q-like quan- 
tity (evolving simply as an exponential) and an 'arrow' 
trajectory always exist. Regarding the system as two op- 
posing teams, we show that the winners always have the 
larger rates-product (e.g., a-c if kakc is larger). When 
the competition is neutral (equal rate-products), there 
will always be a fixed line and two i?-like invariants. 

Following Section II, we now consider the species Xm 
(m — 1, ■ ■ ■ , S) interacting cyclically: 



-t- Xf^ 



(38) 



(with xs+i — xi). The rate equations for the fractions 
Xjn = Nra/N take the form 



(39) 



Of course, J2m-^"^ = 1, so that our configuration space 
is an 5* — 1 simplex (a hyper-tetrahedron). As in Section 
III, a better approach is to cast eqn. (p9l) as 

dt (lnX„) = krnXrn+l - km-lX^^i . (40) 

Despite the non-linearity, we will find it convenient to 
exploit the notation of vectors and matrices. Thus, we 
rewrite this equation as 



where 



5tln_ 



/ fci 
-fci k2 
-fca 





V ks 



KX 



(41) 



-ks 





ks-1 
-ks-i 



(42) 



is an antisymmetric, cyclic matrix and X is a vector with 
elements Xm- Introducing 



£; = 1 1 



1 1 



we see that J2m -^"i = 1 is simply E ■ X = 1. The signif- 
icance of eqn. (|4ip is clear: It highlights the difference 
between having a singular K or not. Such properties are 
best explored by considering systems with odd/even S 
separately. Indeed, the major differences between the 
3SS and 4SS can be understood in the light of this 
even/odd distinction. (Odd-even effects were briefiy dis- 
cussed in fsi!] for systems on a two-dimensional lattice). 



Antisymmetric matrices have a simple property: If k is 
an eigenvalue, then so is —k. Since S is odd, there must 
be at least one eigenvalue which is zero (the rest being 
pairs of opposite signs). Indeed, our K has only one such 
eigenvalue, as shown in Appendix B, so that its null space 
is one-dimensional. The associated right eigenvector, 
satisfies K( — 0. All its elements, (m, can be chosen to 
be positive, since they obey km-iCm-i = kmCm+i- Of 
course, when normalized appropriately, it becomes the 
fixed point population: 



X* 



C 



(43) 



Since IK is antisymmetric, C is also the left eigenvector 
((^•K = 0). Thus, we arrive at dtC-\njl — 0, which leads 
us to define an invariant 

n^l[xt ■ (44) 

m 

Note that ( is determined up to an overall constant, 
which simply corresponds to TZ being fixed (up to an 
overall power). This is the generalization of R in the 
3SS. Together with E ■ X — 1, this invariant defines a 
compact, — 2 dimensional manifold on which the orbit 
must lie. We find it remarkable that ( plays a 'dual' role, 
serving to pinpoint both the fixed point X* and to define 
the invariant TZ. An interesting question naturally arises: 
Does a deeper connection between these apparently un- 
related quantities exist? 

B. Systems with even S 

We may regard such systems as the competition be- 
tween two 'teams,' each with s = S/2 species. In partic- 
ular, this aspect becomes transparent if we reorder our 
species to 

Xi,X3, Xs-l,X2, Xs 

and define 'team variables' 

Y( = X21-1] Zg = X21 



with £ = 1, ...,s. In terms of these, the matrix 
the form 

OR' 



and eqn. (|4T|l becomes 



takes 



(45) 



(46) 
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Of course, the full space and 'team space' have different 
dimensions {S vs. s), but using the same notation (K, M, 
X, Y, etc.) should not cause much confusion. Note that 
eqn. (j46p exposes a simplectic structure that is typical 
in the 'competition of two (sets of) variables.' 

Careful accounting of the s x s matrix M leads to: 











\ 



-ks 





(47) 



'ks-2 ks-i I 



Note that its structure is quite different from K, with 
the (even) odd rates lining the (sub-)diagonal. Whether 
K is singular or not is just controlled by A = det M. A 
straightforward computation leads to 



A= n ^™ n 



odd \ 



regardless of whether s is even or odd. Clearly, A is 
the generalization of A in the 4SS and its sign deter- 
mines which team wins. If A 0, K"-'^ exists so that 



E ■ 



E ■ X = I. Denoting ^ - {R- 



by Um and defining 



exp 



E-K-^hT^] =WX 



(48) 



we see that it evolves trivially: Q (0) e*. Of course, it is 
impossible to form such teams in systems with odd S, 
there is no Q-like variable, and this simple scenario is 
invalid. 

Though obviously related to Q in the 4SS, Q is less 
transparent, since e* obscures the A dependence of our 
system! To find the generalization of Q and to shed light 

on the A — > limit, we exploit the team variables {y , Z 

and consider W, the cofactors of M. Unlike M^-'^, 



(49) 



W = AM" 



is well behaved even when A = 0. Multiplying eqn. (j46 
by 



-W^ 
W 



we have 




= A 



(50) 



(51) 



Now, a tedious computation shows that all the matrix 
elements of W are positive (sums of products of s — 1 
rates, e.g., fcf, + kc in the s = 2 case above). Let us define 
the following positive quantities 



(52) 



which are readily recognized as E - W and E ■ , re- 
spectively. Thus, by taking the scalar product of E with 
eqn. ((?T|) . we are led to the combination 



q-h^-r- h^j = Y[ / Yl (53) 



Q = exp 

which evolves according to 

Q (i) = Q (0) e^* 



(54) 



To help the readers, let us connect these results to the 
4SS case explicitly. There, 




and 



ka kd 

kb kc 



A = A: 



kc kfi 
kb ka 



so that all the conclusions reached in Section III are re- 
covered. 

As in the 4SS, the orbits are, in general, analyti- 
cally inaccessible. However, the 'straight arrow' trajec- 
tory can be determined explicitly, while the time de- 
pendence along it remains the same for any number of 
species. In particular, we begin with the most trivial 
case: s = 1 = S/2: 



dtXi = kiXiX2] dtX2 = -kiXiX2 



(55) 



so that A is just ki. and Q — X1/X2 oc . Since X2 — 
1—Xi, the solution is trivial, being the same form as ([M)) 
with ki in place of w. As t runs to — cxd or -|-oo, it reaches 
the two terminals A"i = 1 or X2 = 1. Turning to general 
S, we seek such a straight line, bridging the terminal 
points {Yi and Zi) in the subspaccs of each team, that 
can serve as a t-dependent trajectory. Making the same 
Ansatz as in eqn. (l30|) . let us consider 

m it) ,Zi[i-h (t)] 

which runs from (Yf,0) to {0,Ze). We wiU also assume 
the same h (t) as in ([M]) . except for uj being replaced by 
f2, a parameter to be determined. Inserting these into 
the left of eqn. (|46p . we find 



E dt In h 
Edtln [I - h] 



E{l-h) 
-Eh 



Meanwhile the right of eqn. (j46|) is 



MZ 
-fA^Y 



E{l-h) 
-fh 



(56) 



(57) 
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where the components of T and :z on the right are, re- 
spectively, 



J2 MenYe and ^ AUeZe 

e I 

Setting expressions (j56l57|) equal, we find 

I 

Exploiting (|49l52p . we have 



(58) 



A 



(59) 
(60) 



Now, all details for our 'arrow' condition are fixed, by 
imposing Yn — 1 etc. The result is 



n = K / E^^' 



(61) 



i.e., the generalization of 

We end with some remarks on the cases with A = 0. 
Although properties of these can be found by taking the 
A limit of the above, let us provide a more direct 
route. Since S is even, a zero eigenvalue of K is doubly 
degenerate and the null space is two-dimensional (at the 
least). Given the form ps)) . we see that this null space 
is determined by the right eigenvectors of M and M"^: 



ie = o = 



We can impose the constraint and define fixed points 
within the space of each 'team': 



Yn* 



/ E ^" ' Zf 




Joining these is a fixed line in the full S —1 simplex, that 
can be parametrized in the same manner as in eqn. (j29p 
of the 4SS: 







7 + 





Z* 



(1-7) 



Meanwhile, the left eigenvectors of K, are expected to 
play a 'dual' role, as in the odd S cases. Of course, 
these can be formed from ff and ^, which are the left 
eigenvector of M and M"^, respectively. The results are 

(^ff, Oj and ^0, ^ , which can now be exploited to define 

the generalizations of / and g, see eq. ([TS]) : rf • In F and 
Not surprisingly, these are intimately related to 



the numerator and the denominator in expression (|53p . 

Note that the sum fj ■ is also an invariant. 

Exponentiating it will produce a quantity which has the 
appearance oiTZ, i.e., products of only positive powers of 
Xm- For systems with all interaction rates being equal, 
such an invariant takes an appealing form: XiX2.-.Xs 
(for any S). Finally, we note that M has no other zero 
eigenvalues (shown in Appendix B), so that there are 
no other invariants. Thus, the orbits lie in a compact, 
5 — 3 dimensional manifold in general. Of course, for 
S" = 4, the 1-dimensional manifolds are simply closed 
loops, for which many properties are explicitly obtained 
in Section III. For S* > 4, we have not found similar 
explicit characteristics. 



VI. SUMMARY AND OUTLOOK 

In an effort to understand better the counter-intuitive 
phenomenon of "survival of the weakest," which was dis- 
covered in recent studies of a well-mixed system involv- 
ing three cyclically competing species, we investigate a 
system with four species. In a previous letter |31| . we 
reported preliminary findings, showing that this behav- 
ior does not persist. Instead, all systems seem to follow 
an intuitively understandable principle: "The prey of the 
prey of the weakest is most likely to go extinct first." In a 
3SS, this species is also the predator of the weakest, and 
so, its early demise leads to the weakest being the sole 
survivor. By contrast, the players in a 4SS form partner- 
pairs, much like in the game of Bridge. As a result, our 
new maxim hints that the weaker pair is not likely to 
survive. All these studies consist of two aspects, each 
of interest on their own. The first is mean field theory 
(MFT), i.e., non-linear dynamics of a set of determinis- 
tic rate equations. The second is the stochastic evolution 
of a finite population of N individuals. For the 4SS, we 
show how the former, our set of rate equations: (|4l)-([71), 
arises from the master equation which defines the latter. 
This article is devoted to an in-depth exploration of the 
MFT, i.e., solutions to the rate equations. 

Since the evolution of the 4SS takes place within a 
tetrahedron, it is more complex than the 3SS case. Nev- 
ertheless, we are able to obtain considerable details an- 
alytically, mainly through the discovery of hidden sym- 
metries. Not surprisingly, far fewer exact results of the 
stochastic aspect are attainable. Instead, most insight 
into this behavior is gained through computer simula- 
tions. Indeed, in some 'extreme' cases, we found that the 
outcomes contradict the MF prediction completely [sij . 
Instead, these can be understood by a complementary 
version of the above maxim, namely, "The prey of the 
prey of the strongest is most likely to survive." A simi- 
lar, in-depth investigation of the second aspect (stochas- 
tic evolution) is beyond the scope here, and those results 
will be reported elsewhere [soj . 

Unlike the 3SS, in which all MFT orbits are closed 
loops (characterized by R) in a plane, the 4SS supports 
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a richer variety of trajectories. The details depend on 
the initial conditions and the rates, especially through 
A — kakc — kbkd- We show that, provided they start 
within our tetrahedron, all orbits spiral, to various de- 
grees, into the a-c edge (A > 0) or the h-d edge (A < 0). 
There is one exception, namely, a special straight line 
bridging the two edges. On it the system moves mono- 
tonically, leading us to refer to it as the 'arrow.' In gen- 
eral, the spirals wrap around this arrow. To describe the 
motion along the spirals and the arrow quantitatively, 
we find a collective variable, Q, which evolves simply as 
an exponential, see eqns. (|16I17I) . Indeed, letting t run 
from — oo to oo, we see that these trajectories all start 
on one edge (a-c or b-d) and end on the opposite one. 
As |A| — i> 0, the spirals become tighter (lower pitched) so 
that, for A = 0, all orbits are closed loops which resemble 
the edge of a saddle. Meanwhile, the 'arrow' becomes 
a line of fixed points, around which all closed orbits en- 
circle. In addition to Q, another invariant emerges. To- 
gether, these two, see eq. completely specify the 
saddle shaped loops. Though the explicit forms are not 
available, many properties of these loops, e.g., maximal 
extent in principle directions, have been found. Finally, 
by extending these considerations to cyclic competition 
of an arbitrary number, S, of species, we can trace the 
origins of these special quantities (invariants, collective 
variables, non-trivial fixed points and lines, etc.) to gen- 
eral properties of an S" x antisymmetric, cyclic matrix, 
K. The elements of K are just the rates, so that all details 
of these special quantities are explicitly known. Indeed, 
a recent study finds that such 'hidden symmetries' are 
present in an even wider class of systems involving com- 
peting species [43 |. Of course, as our scope broadens, 
fewer explicit results are available. 

Let us reiterate that the main forte of MFT is to pro- 
vide reliable estimates of the evolution when the popu- 
lation of every species is large. Relying on continuous 
fractions (and continuous time) , it fares poorly when one 
or more species are near extinction and their numbers 
drop to 0(1). Though it fails to predict all eventual- 
ities of a finite stochastic system, it does give us some 
useful hints. Specifically, the proximity of an MFT tra- 
jectory to an absorbing face should be a good guide for 
the extinction probabilities of the associated species. We 
have not systematically explore this issue and such an 
undertaking seems worthwhile. 

Beyond such immediate questions, we believe there are 
interesting avenues for further research. We end with 
pointing out a few here. The sheets labeled by a, 6, etc. 
are clearly significant, as they correspond to the turning 
points of species a, 6, etc. It seems likely that these are 
good candidates for Poincare sections with which to ana- 
lyze our dynamical system. How often orbits pierce these 
sheets may also serve to classify the spirals by the total 
number of turns, x ('windings'), between the end points. 
Clearly, x is zero for the arrow, while our expectation 
is that this value will be infinitesimal for a spiral in its 
neighborhood. In other words, a perturbative approach 



is likely to be successful. Of course, to be quantitative, 
we should look for a good coordinate system for config- 
uration space, guided by slow and fast dynamics, along 
the lines of action-angle variables in classical mechanics. 
For the neutral cases (A = 0), these have been identified: 
/, 5,6' in eqns. (fT8 | IA7p . In the general case, we already 
have one key variable, Q. It seems quite feasible to find 
two others, perhaps in analog to a radius (distance from 
the arrow) and an angle, 0. Once their i-dependence is 
established, then the classification of orbits we conjec- 
tured would be realized hy x = 4>{^ — oo) 4' — — oo). 
It is even conceivable that a deeper connection with sim- 
ilar concepts in scattering theory exists. 
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Appendix A: Properties of closed orbits 

Quantitative details of a closed orbit can be found 
through its projection onto, say, the pa-pd plane. Us- 
ing A = Aopi^'^'', etc., the (generically transcendental) 
equation for this projection is given by eqn. (|23|) : 



1 



Clearly, this is just one of many equivalent representa- 
tions. We begin with finding the extremal points of this 
closed loop. Exploiting such points, we can transform 
this cumbersome equation into a familiar one, (jA6l) be- 
low. Finally, we will show how the period of these orbits 
can be found, explicitly in certain cases. 

To find the extremal points, let us first consider the 
combination involving pa'. 



F ipa) = AopI/"- + CoPa' 



As F is convex, it has a unique minimum value, F, which 
occurs at 



V 



kaAo 



and corresponds to 

V vl/'^'f' 

A= Ao = 



/V \ V 

In other words, F \Pa\ —F- Similarly, the other combi- 



nation, 



G{pd)^Bo 
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has a minimum, G^ occurring at 
kdBo 



V 

Pd= 



kbDo 



V 

^ D- 



uk^ U-kd rtkd T-^ka 



l/(fca + fcd) 



orbit equation. We begin with the observation that both 

V V 

F {Pa) — F and G {pd) — G he in the range [0, M], where 

A V V V A V 

M =F - F=l- G - F^G - G ■ 



Although neither F nor G are functions bounded from 
above, the constraint F + G = 1 does impose upper 
bounds in this context. Thus, we find the largest value 

A V 

F can take must be given by F^ 1— G • Similarly, we 

A V A A 

have G— 1— F ■ While F and G are unique, each is 
associated with two points, in A and D (or pa and pd). 
Labeling these pairs by A± and D±, we see their signifi- 
cance: They are the extremal points of the orbit, i.e., 

Ae[A_,A+], De[D^,D+] . 

V V 

Let us emphasize that A^,D- y^A,D in general: The 
former is the lower bound for A while the latter is asso- 
ciated with the minimum of F. 

Next, we study the defining equations for the extremal 
points: 

F{A±) = 1-G, G{D±) = 1-F. 

These give rise to eqns. p4ll28p in Section IIIc, repro- 
duced here for convenience: 



A± 
with 



j^^^k./k^ _ 



Jdd: 



-ka/kd 



= 1-Ka 
1 



Tjkd Tjk, 



l/(fea + fcd) 



l/{ka + kb) 



Ja = CoAI'""- 



D 

kd/(ka + kd) 



ka/ika+kb) 



ka/ika 



fe(,/(fca+fcb) 



(Al) 
\a2) 

(A3) 
(A4) 

(A5) 



Let us reiterate: When a species is at its maximum or 
its minimum, the fractions of its opposing partner-pair 
assume the same value. Thus, when A — A±, B takes on 

the same value: [k^K B^'' Dq"" \ while!? 

V 

assumes the value D- 

The other pairs of extremal points {B±,G±) can be 
similarly studied. But, it is much simpler to exploit the 
invariant /, g and to recognize that the extremal points 
are paired, so that 



C± = 



B± = 



Though the solutions to equations like should 
provide some indications of the extent of the saddle 
shaped orbits, our goal is the simple form of the 



This motivates us to introduce the variables 



a — ±i 



F{Pa)'F 



M 



P = ±J G{pd)-G 



M 



V 

where the ± follows the sign of Pm— Pm- Clearly, a,/3 G 
[— 1 , 1] and eqn. (|23|) is transformed into the familiar one 
for a unit circle. 



-/3^^1. 



(A6) 



Undoubtedly, there is a global coordinate transformation 
so that all closed orbits are just circles of various radii. 
But, it may be a challenge to find this transformation 
and we will not pursue further. Instead, let us note the 
possibility of defining an 'angle' 9 by 



9 = tan 



\ 



Gipd)-G 



F{Pa) 



V 

F 



tan 



B + D-G 
A + C-F 



(A7) 



with the understanding that appropriate signs of the 
roots be used for the four phases. It is clear that 9 is 
the only variable in our problem, since we have two in- 
variants and a constraint, see eqns. (|T8l) and ([3]). In this 
sense, since the fractions can be written as (complicated, 
implicit) functions of 9, the evolution is controlled by 
a first order ordinary differential equation dt9 ~ J'{9). 
The solution can be formally written as /p 1 /J- = t and 

so the period, T, around the closed orbit is J^^ l/T. Be- 
yond such formal considerations, we believe that 9 can 
be exploited to formulate our problem with action-angle 
like variables. Such an approach is likely to play an im- 
portant role in the study of the full stochastic model, in 
which 9 serves as the fast variable, while /, g will wander 
slowly due to the noise [29| . 

Let us close this appendix with two remarks. We 
should emphasize that the 'unit circle' is just associated 
with a projection of the saddle-shaped orbit. Though the 
time dependence (i.e., motion around the circle) and T 
remain in general quite implicit, some simplifications do 
emerge if certain pairs of rates are equal, e.g., ka — kd 
(and so, ki, = kc). The period, for example, assumes 
a less prohibitive form than 1/T. Instead of 9, we 
define <j) = \npa and arrive at 



T = 



Ik a 

kb 



4'+ 
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where the end points are associated with A± and are 
given by 



Of course, this class of rates includes the most special 
case, with all fc's being equal, in Q. 



Appendix B: Null space of K and M 

Here, we show that IK has at most two zero eigenvalues, 
so that the dimension of its null space is two or less. For 
odd S, the characteristic equation antisymmetric matrix 
(K) must be of the form 



that this result implies the existence of a single (non- 
trivial) fixed point, 'dual' to a single invariant, see eqns. 
(I43l44| . 

For even S, the rearrangement into 'teams' led us to 
consider M instead of K, see eqn. (|45|) . With no partic- 
ular symmetry, the characteristic equation for M is 

= J]^ (fc2«-l - M) - n 
n— 1 n— 1 

= dctM-ri^ + ... 

with Fi being {kik3...ks-i)J2nKn-i > 0- Thus, the 
null space of M is at most one-dimensional, so that the 
null space of K is either zero- or two-dimensional (for 
even S). 



= det [K - kI] = -CiK + CsK^... - k" 



(Bl) 



where the C's are the co-factors. It is straightforward to 
see that Ci is the sum 



det 



/ 





V 



■det 





\ks 



k2 
fca 


















-fc.ci 



det 






/ 






V 



.. 
.. 

fcl 



-k2 








-ks-1 







-ks-1 
... 

k2 ... 

... 

: 
... 
... 



ks-i 
/ 

-fcs \ 



ks-i 










-ks-2 







ks-2 




which is 



/v2«'4...«'5_i 



.kg 



+ kik^...kg_2 ■ (B2) 
A systematic writing of this expression is 



E 



(S-l)/2 



n 

T! = l 



'^e+2n 



Here, indices greater than 5* are understood to mean 
mod S, of course. Thus, for example, the product in 
the £ = S — I term runs from 5 — 1 + 2 ^ Ito 
S-l + S-l=^S-2, i.e., the last term in (|B2]). Since 
this cofactor is positive, there is just one solution to eqn. 
(jBip with K — 0. Thus, the null space of K is precisely 
one dimensional (for odd S). In the main text, we show 
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